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ABSTRACT 

Earth-based rad ar observations of the rotational dynamics of Mercury 



cih: 

. (iMargot et aLl 120121) combined w ith the determination of its gravity field by 



^ . MESSENGER (ISmith et all l2012f ) give clues on the internal structure of Mer- 

^ cury, in particular its polar moment of inertia C, deduced from the obliquity 

(2.04 ± 0.08) arcmin. 

The dynamics of the obliquity of Mercury is a very-long term motion (a few 
hundreds of kyrs) , based on the regressional m otion of Mercury's orbital a scend- 



ing node. This paper, following the study of iNoyelles fc D'Hoedtl (120121 ). aims 
at first giving initial conditions at any time and for any values of the internal 
structure parameters for numerical simulations, and at using them to estimate 
the influence of usually neglected parameters on the obliquity, like J3, the Love 
number k2 and the secular variations of the orbital elements. We use for that av- 
eraged representations of the orbital and rotational motions of Mercury, suitable 
for long-term studies. 

We flnd that J3 should alter the obliquity by 250 milli-arcsec, the tides by 
100 milli-arcsec, and the secular variations of the orbital elements by 10 milli- 
arcsec over 20 years. The resulting value of C could be at the most changed from 
0.346mi?2 to 0.345mi?2 

Subject headings: Mercury — Celestial Mechanics — Resonances, spin-orbit — 
Rotational dynamics 



1. Introduction 



^The space missions MESSENGER (NASA) and BepiColombo (ESA / JAXA) flBalogh et al, 

20071 ) are opportunities to have a better knowledge of Mercury, in particular its rotational 
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dynamics and its internal structure . The size o f an outer naolten core can be inverted from 
the well-known Peale experiment (jPeald 1 19761 : iPeale et al.l |2002| ) in measuring separately 
the obliquity of Mercury, yielding the polar momentum of inertia C, and the longitudinal 
librations, yielding the polar momentum of the mantle Cm- This last assertion is based on 
the assumption that the longitudinal librations of the mantle are deco upled from the other 
layer s , what can be doubtfu l if Mercury has a significant inner core ( IVeasey fc Dumberry 
2011uVan Hoolst et al.ll2012l ). But this does not affect the determination of C. 



The determination of C from the obliquity is based on the assumptions that Mercury 
is at the Cassini State 1, and that it behaves as a rigid body over the timescales relevant 
for the variations of the obliquity, i.e. a few hundreds of years, due to the regressional 
motion of Mercury's orbital nodes. The C assini State 1 is a dynamical equilibriura in which 
Mercury is in a 3:2 spin-orbit resonance (jPettengill fc Dycelll965l : IColombdll965l ). and the 
small free librations around the equilibrium have been damped with enough efficiency so 
that the obliquity of Mercury can be considered as an equilibrium obliquity. iPeald (120051 ) 
estimated the damping timescale to be of the order of 10^ years for the free librations and the 
free precession of the spin, in considering the tidal dissipation and the core-mantle friction. 
The presence of Mercury near the Cass ini State 1, has been confirmed by Earth-based radar 
observations ( iMargot et al.l 120071 . |2012| ). giving an obliquity of 2.04 ± 0.08 arcminutes. The 



variations of the obliquity cari be considered as small and adiab atic ( iBills &: Comstockl 12005 
PealelboOfil : IBoIs fc R,ambauxll2007l : IP'Hoedt fc Lemaitrel [2008h . 



It is co mmonly accepte d that the obliquity e shou ld be inverted using Peale's formula 

JPealelll969h . that reads as (lYseboodt fc MargotI 120061 ): 



Ctl sin L 



CCl cos t + 2nmW' (|e 



123^3 
16 
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where n is the orbital mean motion of Mercury, m its mass, R its mean radius, e its orbital 
eccentricity, C20 = —J2 and C22 are the 2 most relevant coefficients of the gravity field of 
Mercury, and l and Cl are the inclination and nodal precession rate of the orbit of Mercury 
with respect to a Laplace Plane. The Laplace Plane is a reference plane, determined so 
as to minimize the variation of the inclination l. In this way, l and ^ are assumed to be 
constant, and also e as a c o nsequence. There are s e veral ways to define the Laplace Plane 
jYseboodt fc MargotI [2OO6I : IBoIs fc Rambauxl I2OO7I : IP'Hoedt et aDl2009h . and so l and ^ 
depend on the chosen definition. 



In order to bypass the uncertainty on the Laplace Plane, iNoyelles fc D'Hoedtl (120121 ) 
proposed a Laplace Plane-free study of the obliquity, using averaged equations of the rota- 
tional motion, averaged variations of the eccentricity, inclination and the associated angles. 
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and a quasi-periodic representation of these quantities. This yields a quasi-periodic repre- 
sentation of the equihbrium obhquity that is very close to th e Cassini State over the domain 
of vahdity of the JPL DE 406 ephemerides (IStandishlll998l ). i.e. 6,000 years between JED 
0625360.50 (-3000 February 23) to 2816912.50 (+3000 May 06). The inertial reference frame 
is the ecliptic at J2000. This approach allows to consider the small variations of the orbital 
elements. 

The spacecraft MESSENGER c urrently orbiting a round Mercury has recently given us 
accurate gravity coefficients (TabU]) (ISmith et al.ll2012l ). this is the opportunity to refine our 
model of the obliquity and to consider usually neglected effects hke the tides or higher order 
gravity coefficients. 



Table 1: The gra vity field coefficients of Mercury derived from MESSENGER data 
( ISmith et al.ll2012l ). These are unnormalized coefficients while Smith et al. give them nor- 
malized. 



C20 = 


-J2 


(-5.031 ±0.02) X 10-5 


C21 




(-5.99 ±6.5) X 10-^ 


S21 




(1.74 ±6.5) X 10-s 


C22 




(8.088 ±0.065) X 10"^ 


S22 




(3.22 ±6.5) X 10-^ 


C30 = 


-J3 


(-1.188 ±0.08) X 10-5 


C40 = 


-J4 


(-1.95 ±0.24) X 10-5 



We first present how we average the equations ruling the rota tional dynamics of Mercur y 
(Sect 12]). This averaging is more accurate than the one given in (INoyelles fc D'Hoedtll2012l ). 
i.e. it is done to higher order in eccentricity. Then we explain how we derive 2 new formulae 
for the obliquity, one analytically (SectJH]), with an approach slightly different than Peale's, 
and one numerically (SectJH). Then the reliability of these formulae are tested and additional 
effects are discussed (SectJS]). The influence on the interpretation of the internal structure 
of Mercury is flnally addressed (SectE]). 



2. The dynamical model 



The basic model behind the spin-orbit dynamics of the Sun- Mercury system is given in terms 
of the Hamiltonian: 
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H = Hk + Ht- {G,Mm)VG, (2) 

where Gc is the gravitational constant, M is the mass of the Sun and m ~ 1.6 ■ 10^^[M] is 
the mass of Mercury. The symbol Hk labels the unperturbed Kepler problem: 

Hk - , (3) 

with the constant /x given by /i = Gc{M + m), and Lq = rriy/JId is the orbital angular 
momentum, which depends on the semi-major axis of Mercury a ~ 0.387[Af/]. Ht defines 
the free rotational motion of Mercury: 

where A < B < C are the principal moments of inertia, G is the norm of the angular 
momentum G, L = Gcos{J) is the projection of G onto the polar figure axis of Mercury, 
with the angle J usually called the wobble, and I is the conjugated angle to the action L. 



2.1. General form of the potential 



The gravitational interaction of the orbital and rotational dynamics is given by the potential 
Vg. I t can be expanded into sp herical harmonics, and takes the form, after (j Cunningham 
1970l : iBertotti fc Farinellalll99ol ): 



G 



EE 

n=0 m=0 



P„m(sin ip) {Gnm cos(mA) + Sr,m sin(mA)) 



(5) 



where R ~ 2439.7[fcm] (lArchinal et al.ll201ll ) is the radius of Mercury, r is the distance 
of the center of mass of Mercury from the center of mass of the Sun, = -Pnm(w) are 
the associated Legendre polynomials, which are defined in terms of the standard Legendre 
polynomials P„ = Pn{u) by: 



Pnm{u) = (1 - n") 



2\m/2 dPn(u) 



The angles (v?. A) are latitude and longitude with ip G (—90°, 90°), A G (0,360°), and the 
Cnmi Snm are the Stokes coefficients, with m < n and n, m G N. The standard convention 
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is to define Coo = 1 and Sno = 0. The notation J„ = — C„o is also used for the remaining 
zonal terms, with m = 0. In addition, by the proper choice of the coordinate system through 
the center of mass of Mercury, the first order coefficients Cio, Cu and vanish, the same 
is true for C21 = S21 = if the axes of figure are aligned with the main axes of inertia. 
The effect of the perturbation on the rotation is therefore proportional to the size of the 
remaining zonal (m = 0), tesseral (m < n) and sectorial (m = n) terms. The remaining 
second order coefficients, C20 = —J2 and C22, can be related to the principal moments of 
inertia by: 



Jami?^ = C 



(A + B) 



C22mR^ 



B-A 



It is a common practize to introduce the normalized polar moment of inertia c through the 
additional equation C = cmR^. 

The contributions Hk-, Hpt and Vq are given in different reference frames. To match them 
we aim to express both in an inertial frame, say e^. For this reason we introduce the unit 
vector (x, y, z), pointing from the center of mass of Mercury to the one of the Sun, which is 
defined in the body frame, say 63, in terms of the angles {ip, A) by the relations: 



X = cos ip cos A, y = cos (p sin A, z = sin if 



(6) 



Together with the definition of P„ in terms of the sum 



l)>'{2n-2k)\ , 

ill 



2" ^ k\{n - k)\{n - 2k)\ 



we find for u = simf the explicit form for Pnm, which are given by: 



[(n-m)/2] 



Pnm{sm^) = ^ COS"^ {if ) ^ 



fc=0 



{-lf{2n-2k)\ 
k\{n — k)\{n — 2k — m)\ 



Together with the formulae by Vieta (see e.g. iHazewinkell ( 1200 ll )) 



sm 



n— 2fc— m 



(7) 



sm 
cos 



(mA) = J2 



k=0 



^ k ^ ■ m-k A sin / 1 , , , 

COS Asm™ "A {-{m-k)'jT 
k / COS \ 2 



we find from Eq.([5]), Eq.(l6]) and Eq.(l7]): 
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1 °° " 

For n = 2 the terms turn out to be: 

^20 = ^ (3^' - 1) , £21 = Sxz , (2:22 = 3 {x^ - f) , 
and for &2m. we find: 

©20 = , ©21 = 3p , ©22 = 6xy . 
The third order contributions, with n = 3 are 

^30 = Iz {5z^ - 3) , C31 = {5z^ - 1) , 5:32 = I5z {x^ - f) , C33 = 15x (x^ - 3^^) 
for (Esrn, and 

©30 = , ©31 = (5^2 _ 1) ^ = sOxyz , ©33 = -15 {f - 3a;2) , 

for ©3m- In our study we are also interested in the fourth order contributions. For n = 4 we 
find for £4^ 

£40 = ^ (3 - 30^2 + 355^) , 2:41 = hz {7z^ - 3) , €42 = {z^ + 2f - l) {7z^ - l) , 

o Z Z 

£43 = I05xz {x^ - 3f) , £44 = 105 {x^ - 6x^f + f) , 
and the terms ©4^ turn out to be of the form: 

©40 = 0, ©41 = (75' - 1) , ©42 = 15£y (752 - 1) , 

©43 = -105y - 3x2) ^e^^ = 420xy{x^-f) . 

Note, that in the above expressions, (tnm and &nm, it is possible to express one of the 
[x, y, z) through the two others of them by making use of the relation Eq.([6]), i.e. the fact 
that x^ + y2 + = 1. 
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2.2. Matching the reference frames 

To express the unit vector {x,y,z) in cq we make use of the usual Andoyer angles {l,g,h) 
together with the angles {J,K); the angle J was already defined above, the angle K enters 
the projection of G on the inertial z-axis by H = Gcos{K). 

Let us denote by cq the inertial, by ei the orbital, by 62 the spin and by 63 the figure frame 
of references, respectively. The unit vector (£, y, z), given in 63, can be expressed in the spin 
frame 62 by 



{x,y,z)^^ = Ml ■ {x,y,z)^^ 

with the matrix Mi given by: 

Ml — diag (cos(g') cos(Z), cos(g') cos( J) cos(Z), cos( J)) . 
The transformed vector can be expressed in an inertial frame eo using the transformation 

y> ^)eo = ^2 • {X, y, Z)^^ , 

with 

M2 = diag (cos(/i), cos(/i) cos(X), cos(i^)) . 

On the other hand, if the position of Mercury is given in the orbital frame ei it is also given 
in the inertial frame eo through: 

(a;,^,i),^ = M3-(cos(/),sin(/),0) , 

where / is the true anomaly of Mercury and the matrix M3 is given by: 

cos(w) cos(f2) - cos(z) sin(u;) sin(f2) - cos(i) cos(u;) sin(n) - cos(n) sin(a;) sin(z) sin(f2) 
-^3 = I cos(i) cos(f2) sin(a;) + cos(w) sin(f2) cos(i) cos(a;) cos(f2) — sin(w) sin(f2) — cos(f2) sin(i) 
sin(i) sin(a;) cos(w) sin(i) cos(i) 



Here a;, ^, i denote the argument of perihelion, longitude of ascending node, and inclination 
of Mercury, respectively. Prom the decompositions 
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M,=R^{g)Ri{J)R3{l), M2 = R3{h)Ri{K), M3 = /?3(^2)i?i(0^3M 
we are able to find 

{^,y,z)^^ = Rm ■ (cos(/),sm(/),0) , 
where we introduced the rotation matrix Rm being of the form: 

Rm = R3{-l)Ri{-J)R^{-K)R^{-h)R^{a)Ri{i)R^{i^) ■ 
Together with the relations 



a 
r 



00 



1 + 2^ J^(z/e) cos (z/M) , 



u=l 



cos(/) = 2 ~^ J^{ve) cos(z/M) - e 

^ u=l 



sm 



(/) = 2^r^^ 



^ dJ„{iye) sin(z/M) 



de V 

u=l 



where are the Bessel functions of the first kind we are able to express the potential Eq.dH]) 
in terms of the rotational and orbital elements only: 

V = V{l,g,h,J,K,a,e,i,u,M). (9) 

It can also be expressed in terms of suitable action angle variables {k, Li), with i = 1, . . . ,6, 
using the set of modified Andoyer variables 



(/i, k, h) = {l + g + h, -I, -h) , 
(Li, L2, L3) = {G, G-L, G-H) = (G, G{1 - cos( J)), G{1 - cos{K))) (10) 
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for the rotational motion, and by making use of the classical Delaunay variables 

(/4, k, k) = (M, 00, a) , 

(L4, L5, Lq) ^ (jjo, L4VI - e^, L5Cos(i)^ (11) 

for the orbital dynamics. In this setting the potential can be written in the form V = V{1, L) 
with / = (/i, . . . ,Iq) and L = (Li, . . . , Lq). 

2.3. Simplifications and assumptions 

The rotation period of Mercury, about = 58.6[d], and the orbital period around the sun. 
To ^ 87.9 [d] lie close to the 3 : 2 resonance (2To ^ 3T,.). It is thus desirable to find a much 
simpler dynamical model, which reproduces the qualitative dynamics close to the resonance. 
We first introduce the change of coordinates 

with a — ((7i, . . . , (Je), S = (Ei, . . . , Eg) defined by the generating function 5'3:2 of the second 
kind 

'S'3:2 = El ^/i — -I4 — I5 — Iq^ + E2Z2 + E3 {I3 + Iq) + E„Z«; . 

In this setting the relevant resonant dynamics can be easily described in terms of the variables 

3 _ _ 

Ci — h — -I4 — h — h , cr2 — I2 , crs — I3 + Iq , 

El = Li , E2 = iv2 , E3 = L3 , 
while the remaining variables become: 

2E4 = 2L4 + 3Ei , E5 = L5 + El , Eg = Lq + El — E3 , 

and (Tj = li with i — A,5,6. In our first approach we aim to construct a simple resonant 
model, valid only close to exact resonance, by making use of the following assumptions: 
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i) we neglect the wobble motion of Mercury; we therefore set J = 0, which implies L2 = 
and reduces Eq.([3]) to: 

(12) 

ii) we neglect the effect of the rotation on the orbital dynamics, i.e. we investigate the 
dynamics on the reduced phase space 

dai _ dH dT.i _ dH 
dt dT^i ' dt dai 

with z = 1,3, therefore assume that the orbital parameters a,e,i,u,^,M are known 
quantities, and Sj = Sj(t),crj = ai{t) with z = 4, 5, 6 act as external time- dependent 
parameters on the dynamics of the reduced phase space. 

iii) we neglect short periodic effect^ and replace the potential Vq by its average over the 
mean anomaly of Mercury M = I4 = a^. 

27r 



{V)^{VG)^^ = ^JJvG{cT,J:)da,. 



As a result the averaged potential {V) becomes independent of a^, or the mean anomaly 
M, and thus the conjugated action S4 becomes a constant of motion, and as a conse- 
quence the semi-major axis a is a constant too. 

iv) we assume the presence of a perturbation leading to an additional precession of the 
nodes with constant precession rates, say w, f2 7^ 0. In this setting we find for the 
remaining phase state variables, connected to the orbital motion, and ctq- 

da5 . dcTQ ■ 

Note that since as = 0-5 (t) = ut + uq and = a^it) = Clt + ^0 the time dependent 
generating function S'3:2 leads to 

at 



^ The effects within time scales, which are smaller than the revolution period of Mercury around the Sun. 
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which we have to add to our new Hamiltonian. The final resonant model takes the form 

Hj = ho + {V) , (13) 

with the new fundamental part: 

ho = Hk + Hr-Eilo + {Es-E,)^ (14) 
and where Hk transforms into the new expression: 

Hx = n ■ 

(2S4-3Si)^ 

We derive from Eq.f fMl) the unperturbed angular frequencies: 



"^ = 9ST = C -2(S4-is03-"-^ = ^^-2"-"-^' 
dho . 

since 

Here n is the mean motion of Mercury. For cj = ^ = the spin-orbit resonance of Mercury 
translates into the commensurability: 

2gi = 2/1 — 30-4 ~ , (Ts — CTq ^ , 

while for ^ 7^ small frequency corrections due to the potential (y) have to be taken 
into account. 



2.4. Workout of the time dependent resonant model 

In the following discussion we limit our investigation to the contributions of the averaged 
potential, in which the Stokes coefficient in Table [T]are bigger than the threshold 10~^, which 
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turn out to be C20, (^22, C^q, C40. According to this simplification the averaged potential (V) 
can be split into the form: 

{V) = C20 {V20) + C22 {V22) + C30 {V30) + C40 {V40) , 

where we used the notation {Vnm) to indicate the terms proportional to the coefficient Cnm- 
For the ongoing investigation we also need to separate the individual terms into purely 
resonant terms, just depending on (Ji^cr-^, and time dependent resonant terms through the 
presence of the additional angles = I5 = U!,aQ = Iq = with cu = uj{t),^ = ^{t). We 
use the short-hand notation: 

{Vnm) = {V„m) {<^1, f^S, h, k) = {v-am) (o^l, f^^) + («nm) ((^l, (^S, h, k) (15) 

to indicate the dependency on the angles of the various terms. 

The term (1^20) = 

The expression proportional to C20 in the unaveraged potential Vq takes the form: 

which becomes after the average over the mean anomaly, and expanded up to 4th order in 
the orbital eccentricity e: 

= _:5!li±^±i^ + [2U cos (^3) + [3]2o cos (2<J3)) + 0{e') . 

The coefficients \j]2o, j — 1,2,3 are depending on Ei, E3 through the inertial obliquity 
K, as well as on the orbital inclination i, and can be found in the Appendix. Note, that 
(V20) = (^^20) (cs) does not depend on ai and is free of the angles ^5,^6 up to 0{e)^ (i.e. 

(«20) = 0). 



The term (^22): 
The expression proportional to C22 in Vq is of the form: 
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The averaged term up to O(e^) becomes: 

{V22) = ((^22) ((71, (Jz) + ^22) ((71, (73, ^5)) + O(e^) , 

with 

{V22) - ([1]22 COS (2ai) + [2]22 cos (2cti + ag) + 

[3)22 cos (2(7i + 20-3) + [4] 22 cos (2(7i + 3(73) + [5] 22 COS (2(7i + 4(73)) , 

where we collect the terms depending on the orbital eccentricity e by: 

Ei = e (-56 + 1236^) . 

The time dependent part of the expression is 

(1^22) =318e^([6] 22 cos (2(7i + 2/5) + [7] 22 COS (2(Ti + (73 + 21^) + 

[8]22 COS (2(Ti + 2(73 + 2/5) + [9]22 COS (2(7i + 3(73 + 2/5) + [10]22 COS (2(Ti + 4(73 + 2/5)) , 

where the coefficients [1]22 ■ ■ ■ [10]22, depending on the resonant actions, are given in the 
Appendix. Note, that (^22) enters a term with the only resonant argument 2(7i, while the 
other Fourier modes are of the form 

2(7i + A;cT3, A; = 1,2,3,... 

up to order 0{e*). The terms (^22) = ('"22) (73, l^) do not depend on Iq and contain terms 
of the form 

2(71 + A;(73 + 2^5, A; = 1,2,3,... 
with the additional argument 21^. 



- 14 - 



The term (1^30) : 
The third order term proportional to C30 turns out to be: 

R"^ 1 R"^ / o \ 

After the averaging no pure resonant terms survive ((t'ao) = 0) and (V30) takes the form 

^3 

with 

(uso) = E2 ([Ijao sin (k - SfJa) + [2)30 sin (^5 - 2(73) + 
[3]3o sin (Z5 - (73) + [4]3o sin + [5)30 sin (Z5 + (T3) + 
[6]3o sin (k + 2(73) + [7]3o sin (k + 3(73)) ■ 

In the above expression the [7] 30, with j — 1, . . .7, label the third order coefficients, depending 
on the resonant actions, which are collected together in the Appendix. Moreover, we define 

E2^e{2 + 5e') 

to collect the contributions, which depend on the eccentricity. Note, that the time dependent 
terms (u^o) are again independent of Iq (they only depend implicitly on it through the 
definition of (T3). We also conclude, that C30 does not contribute to the time independent 
resonant model, at the first order of masses approximation. 

The term {V^q): 

The expression proportonal to C40 of {Vq)^ originates from the term 




and can again be split into 
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(^4o) = Y024^ ^^^^ ^ ^^^^ ■ 

Here the time-independent part is given by the expression: 

{V40) = 4£;3 ([l]4o + [2]4o COS ((73) + [3]40 COS (2(73) + [4]40 COS (30-3) + [5]40 COS (4(73)) , 

where 

£;3 - (1 + 5e2) , 
and the time dependent contributions are: 

{U40) = E4 ([6]40 COS (2I5) + [7]40 COS (2^5 - 4(73) + [8)40 COS (2^5 - 8(73) + 
[9]40 COS (2^5 - 2(73) + [10]40 COS (2/5 - (73) + [11]40 COS (2^5 + (73) + 

[12]40 COS (2/5 + 2(73) + [13]40 COS (2/5 + 3(73) + [14]40 COS (2/5 + 4(73)) , 

where 

£;4 = (2 + 7e^) 

(the coefficients [. . .]4o can be found again in the Appendix). Note, that (^40) just depends 
on (73, while (^40) depends on (73, 2/5 but not on Iq. 

To summarize, we find: 

(i) The effect proportional to C20 (or J2) is of order R^/a?, is time independent, and 
depending on the resonant angle (73 only. 

(ii) The effect proportional to C22 can be split into time dependent as well as time inde- 
pendent terms: the time independent terms are proportional to ei?^/a^, while the time 
dependent contributions are of order e'l^jc?. The former contains a Fourier term 
depending just on the resonant argument 2(7i, while the later is depending on integer 
combinations of (7i, (73 and Z5 only. 



(iii) There is no time independent effect proportional to C30 (or J3) on the averaged dy- 
namics. Only time dependent terms, which are of order ei?^/a^, and depending on 
integer combinations of (73 and Z5 contribute to it. 
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(iv) There is an important, time independent contribution, of order i?^/a^, which is pro- 
portional to C40 and just depending on the resonant angle a^. The time dependent 
part of the potential is of order e^R^/a^ only. 

The preceeding list shows, that the time dependent effects are smaller than the time inde- 
pendent. For the numerical integration of the averaged system (see SecJl]) we are going to 
use the full potential of the form Eq.f lT^ . while for the analytical study we are going to use 
the time independent part of the potential only. 



In this section, to obtain a simple analytical formula, we neglect the resonant terms in {V), 
which depend on time through ^5 = u{t), Iq = Q{t), and investigate the long-term dynamics 
close to (Ji = (Js = 0. The potential {V) = {V) (ui, 0-3) reduces to 



3. Analytical treatment of the Cassini State 1 



{V) (cTi, 0-3) = C20 {v2o) (c^a) + C22 {V22) (o"i, 0-3) + Cio {vio) (0-3) , 



the integrable parts reduce t< 



and the new Hamiltonian model becomes: 



Hio = Hk + Hr+ (S3 - Si) ^ + C20 {V20) + C22 {V22) + C40 {V40) 



(16) 



The requirement to remain at the equilibrium is that Si = S3 
set of equations: 



and translates into the 




(17) 



^ There is no need to keep u in the ongoing calculations. 
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The system can be solved for Si = Si^,, and S3 = S3*, which imphes the 'equihbrium norm' 
of the angular momentum C and the 'equilibrium obliquity' i^^,, which comes from the 
equations Si^, = G^, S3* = Si* (1 — cos {K^)). The physical interpretation of the equilibrium 
solution is the following: ai = means that the axis of smallest inertia points (on average) 
towards the Sun, (T3 = ensures that the node of the equator of Mercury is locked with the 
node of its orbit. While G* defines a small correction of the unperturbed spin frequency, 
the angle defines a specific value of the inertial obliquity. It translates to the usual called 
obhquity e by the relation 



cos(e) = cos(i) cos(i^') + sin(i) sin(i^') cos (a3 — a^) , (18) 

which reduces to 
for 0-3 = cTg = 0. 

The contributions (fnm) depend on (Si, S3) through Sk = s/^ (Si, S3) , cj^ = c/^(Si,S3), 
since from Si = Li, S3 = L3 we find Si = G, S3 = Si(l — cos(K)). We aim to express 
Eq. llTTl) in terms of (G, K) and thus have also to express the derivatives d /9Si , d /5S3 in 
terms of (G, K) too. A simple calculation shows for sk and ck and its derivatives: 

dcK ^ — Ck dcK 1 dsK ck — 1 Osk 1 



aSi G ' 9S3 G ' aSi GtK ' 9S3 Gt 



K 



where we have introduced tx = tan(i^), and used the relation sk^ + ck^ = 1 and < K < 
7r/2 to simplify the expressions. A long but straightforward calculation shows, that Eq.( lT7|) 
can be written as: 

/l(G, K) = + + I — ([IJZ/Cso + [Af^22j + ^[3J/C40 ) = ' 

/2(G, K) = ^ + (^^ ([4];G2o + [5]/G22) + ^[6]/C^4o) = , (19) 



with 



[1]/ = (8 + 12e' + 15e^) S2i^-K)tK/2 , 
3 

[2]/ = -gC (-56 + 1236^) cl^K)/2Hi'K)/2tK/2 , 

[3]/ = (8 + 40e2 + 105e^) (2s2(i-x) + 7s,^,_k)) 
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and 

[4]/ = ^ (8 + 12e2 + 15e^) S2(^^-k)Isk , 
3 

[5]/ = — e (-56 + 1236^) {2cKSi + S2(i_i^)) /sj^ , 
15 

[6]/ = (8 + 40e2 + 105e^) (2s2(._i^) + 7s4(.-i^)) /s;^ . 

Note, that we can use the second of Eg. (1191) to ehminate G from the remaining equation to 
get: 

F{K) = - (y + tlcK^ + GcMm ([IJ^Cso + [2]fC22) + ^[3]fC4o] = , (20) 



with 



= -7j7^ (8 + 12e2 + 15e^) S2i^-K) , 
[2]f = -777777-^ ;3e (-56 + 1236^) , 



QAC^SKS{i-K)/ 



2 



15 

Since the observed obhquities are usually rather small we substitute Eq.f lTS]) for = = Q 
in Eq. (l20l) and expand around e = up to first order. Moreover we make use of the relations 
C = cmR^ and n^a^ ~ CM to get rid of some constants. Within these approximations we 
arrive at the simple formula: 



1 + -— cos(z^ 
3 n 



(21) 



cCl sin(i) 



n 



2C22 He - ^e3) - C20 (1 



3p2 
2^ 



5 I 25 2 
2^2 ' 16 



I ( ^ ) csin(z)2 



Notice, that here the parameter ^ is assumed to be negative, and that i is defined with 
respect to the plane, which minimizes the variations of the inclination. The formula coincides 
for C40 = with Peale's formula up to 0{^/n)'^ with the exception of the term c^cos{i) in 
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the denominator ( Yseboodt fc Margot 20061 ) However, a comparison of Eq.f l2ip with two 



different versions of Peale's formula clearly show the validity of the approach (see Fig{T]): 
within the interval c G (0.3, 0.4) the predicted equilibrium obliq uity at Cassini state 1 frora 
Eq. fl^Tl) lies well within the value predicted with the formula in JVseboodt fc MargotI [20061 ) 



and ( lPealelll98ll ) with a difference of about lOe^ ~ 600mas. Note that the parameter C30 
is absent in this simple averaged model (there is however a time dependent effect as we will 
see below), the influence of C40 on the denominator is of the order of 0{R/af' ~ 4. ■ 10~^, 
and 0{^/n)^ ^ 8. • 10^'' does not modify the results. For the parameters given in Tabled] 
the equilibrium obliquity for Mercury turns out to be e = 2.02'. 

We conclude the section with a short summary of the assumptions, which were made to 
obtain Eq. fl21l) : i) zero wobble, ii) only the coupling of the spin on the orbit was taken 
into account, iii) short periodic effects are neglected (average over mean orbital motion and 
time), iv) the orbital parameters are kept constant but then we assume a constant rate of 
regression of the ascending node. The resulting accuracy of the formula can be quantified 
as follows: 1) for the average a 4th order expansion in eccentricty e was taken into account 
(error O(e^) ^ 3. ■ 10""^), 2) time dependent resonant terms are neglected (with this we set 
de / dt = 0) , 3) a. first order expansion in obliquity e allowed to make the formula explicit and 
induced an aditional error of O(e^) ^ 1. ■ 10"'^. 



This formula and Peale's formula as well are lacking of the fact, that the obliquity changes 
with time not only because f2 is not a constant but because the 'equilibrium' conditions 
c"i = = cannot be fulfilled for all times. As we have seen, even for constant precession 
rates, the Hamiltonian is not a conserved quantity, i.e. for a Hamiltonian of the form 
H = H{T,i, S3, cTi, (T3, t) we cannot deduce that cxi = (72 = for all times, since 

9i/(Ei,E3,(Ti,(T3,t) 

and therefore a\^2 are subject to change too. A second (and probably more visible) concern 
of formulas of the type Eq. fl?!]) is the fact, that the resulting e strongly depends on the choice 
of the reference plane to which the orbital inclination i and in which the average precession 
rate ^ are defined. To optimize the result, the reference plane should be the plane to which 



•^For small {tl/n) it is possible to simplify the calculations by setting G — UgC ~ ^nC in the second of 
Eq. (|19p . expanding up to 1st order in e close to e = 0. The resulting formula for e coincides with Eq. (|2ip 
up to 0{fl/n)'^. The additional terms are therefore stemming from the small spin frequency correction, not 
taken into account in (Peale 1981). 
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the variations in i become minimal, which can be seen as a generahzation to the so-called 
Laplace plane (not to be confused with the invariant Laplace plane, the plane normal to 
the angular momentum of the complete system). In the next Section we present a possible 
solution to these kind of problems. 



The secular variations of the orbital elements 



We use a similar approach as in (iNoyelles &: D'Hoedtll2012l ). In a re ference frame based 
on th e ecliptic at J2000, we define averaged eccentricities and inclinations ( iNoyelles fc D'Hoedt 
2012h as: 



h{t) = -7.76651 X IQ-^H^ + 1.43999 x lO^^t + 0.200722, (22) 

k{t) = -2.31417 X lO^^V - 5.52628 X 10"^t + 0.0446629, (23) 

p{t) = 2.38036 X lO^^^t^ - 9.03918 x IQ-^H^ - 1.27635 x IQ-H + 0.0456355, (24) 

p{t) = -1.04673 X IQ-^H^ - 1.27792 x lO^^t + 0.0456362, (25) 

q{t) = 2.52407 x IQ-^H^ - 1.06586 x 10"^^^^ + 6.54322 x 10~h + 0.0406156, (26) 

q{t) = -1.21729 X lO^^H^ + 6.52656 X 10"^t + 0.0406163, (27) 

with h{t) = e(t) smzu{t), k{t) = e(t) cos zu(t), p{t) = sin (^^^ sin^()f:) and g(t) = sin (^^^ cos^(t), 
the time origin begin J2000 and the time unit being the year. / is the orbital inclination of 
Mercury defined with respect to the ecliptic at J2000. These are fits over the duration of 
the JPL DE406 ephemerides, i.e. 6,000 years. For the inclination variables, 2 fits are given: 
the third-order fit (EqJ2Hand Eq J2Bl) is more accurate, but the second-order one is easier to 
extrapolate into a trigonometric decomposition. 

From these formulae we extract trigonometric expressions: 



h{t) = 0.1990903983 sin ti7i(t) +0.01094807206 sin ti72(t), (28) 

k{t) = 0.1990903983 cos ti7i(t) + 0.01094807206 cos ti72(t), (29) 

p{t) = 0.06094690052 sin (t) + 0.01442538649 sin fi2(t), (30) 

q{t) = 0.06094690052 cos fii(t) + 0.01442538649 cos 1^2 (t), (31) 

with 



= 2.852011398 x lO'^t + 1.30845314198, 



(32) 
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tU2{t) = 4.767836272 x IQ-H + 2.26085090227, 
fii(t) = -2.298222197 X 10~^t + 0.60658814513, 
fi2(t) = 1.340719884 X 10-^t + 2.28580288184. 



(33) 
(34) 
(35) 



Using trigonometric series make the orbital solutions easy to extrapolate without di- 
vergence. We extrapolate them over several millions of years to optimize their numerical 
identification, and we assume that the solution is close to the real one over the validity of 
the DE406 ephemerides. This assumption will be checked a posteriori (Sect l5.1"j) . 

We can now plug these new orbital elements in the e quations of the ave raged rotational 
dynamics of Mercury. Using Laskar's frequency analysis (jLaskarlll999l . l2005l ). we express the 
ecliptic obliquity K and the resonant argument (Xs with a quasiperiodic decomposition such 



as 



K{t) = /(t) + ^ajcoswit + 

(73 (t) = bj COS COjt + (f)j, 



(36) 
(37) 



Qi and bj being real amplitudes, Uij frequencies, and (pij phases at t = 0. The frequencies 
can either come from the forced motion, and so are combinations of the ones present in the 
orbital motion (Eqj32]to Eql35l). or are due to free librations, that are expected to be damped. 
Their presence in the numerical outputs comes from a non optimal choice of the initial 
conditions. Our first run with C20 = -5.031 x 10-^ C22 = 8 x 10^^ C30 = -1.188 x 10~^ 
C40 = -1.95 X 10-^ and C = O.SbmR'^ yields the TabU 



The frequency analysis giving this table has been performed over 7.95 Myr with 4,096 
points equally spaced by 1,942.5 years. The fr ee librations should have pe r iods of ^ 15 years 
in lon gitude and ~ 1,000 years in obliquity fjP'Hoedt &: Lemaitrd |2004J : iRambaux fc Bois 



2OO4J ). these periods should appear aliased, while the forced perturbations should not since 
their periods should be bigger than 10,000 years. We the n optimize the initia l conditions 
iteratively in removing the free librations, as described in ( Noyelles et al.l 2013 ). These free 
oscillations act as a noise in the determination of the forced ones, that explains for instance 
small discrepancies in the periods. Refining these initial conditions improves the accuracy 
of the determination of the forced oscillations. 
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Table 2: An example of frequency analysis of the numerical outputs of our system. 



N 


Amplitude 


Period 


Identification 




(arcmin) 


(y) 




K -I 


1 


1.9755 


oo 


<K-I> 


2 


0.2682 


172,665.2 




3 


0.0592 


86,332.6 




4 


0.0254 


264,530.5 




5 


0.0132 


60,999.1 


2uJi - 2VLi 


6 


0.0121 


57,555.3 


3^^2 — 31^]^ 


7 


0.0026 


4,167.1 


free 


8 


0.0025 


43,166.2 


4^2 - 4fii 


9 


0.0025 


4,302.6 


free 




1 


6.2509 


172,665.2 


1^2 — ^1 


2 


1.4683 


86,332.6 


2^^2 — 2^^]^ 


3 


0.3462 


57,555.1 


3^^2 — 31^]^ 


4 


0.1112 


60,999.0 


2uji - 2Vti 


5 


0.0816 


43,166.3 


AVL2 - AVLi 


6 


0.0400 


497,291.8 


W2 — '^l + ^2 ~ ^1 


7 


0.0396 


104,473.0 


Wi — W2 + ^2 ~ ^1 


8 


0.0261 


45,075.0 


2wi — 3r2i + 


9 


0.0214 


4,167.1 


free 


10 


0.0209 


4,302.6 


free 



4.1. Fitting the initial conditions 

The goal is to find relevant initial conditions for any set of interior parameters (C205 ^^22, C'so, ^^405 C) 
consistent with the observations and the theory. For that, we consider 56 different cases 
(TabJa}. 
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Table 3. Our 56 sets of possible interior parameters. 



N 


Gn (xlO-^) 


C22 (xlQ-^^) 


r,o(xio-^^) do 


(xlQ-'^) 


C/(rriR^) 


n 
u 


(Til 


o.Uoo 


1 1 aa 

— i.iOO 


— i.yo 


n Qc; 
u.oo 


i 


— O.Uoi 



0. 


1 1 aa 

— i.ioo 


— i.yo 


n Qc; 
u.oo 





— O.Uoi 


a no 
o.uz 


1 1 aa 
— i.ioo 


— i.yo 


u.oo 


Q 
O 


— O.Uoi 


a n/i 

O.U4 


1 1 aa 
— i.ioo 


— i.yo 


n Qc; 
u.oo 


A 


— O.UOi 


a nfi 

O.UD 


1 1 aa 

— i.ioo 


— i.yo 


u.oo 





— O.UOi 


8 na 

O.Uo 


1 1 aa 

— i. ioo 


— i.yo 


u.oo 





— O.UOi 


a 1 
0. i 


1 1 aa 

— i. loo 


— i.yo 


u.oo 


7 


— O.UOi 


a 1 
0. iz 


1 1 aa 

— 1. ioo 


— i.yo 


u.oo 


o 

o 


— O.UOi 


a 1 /I 

0. i4 


1 1 aa 

— 1. loo 


— i.yo 


u.oo 


Q 

y 


— O.Uoi 


a 1 
0. iu 


1 1 aa 

— i. ioo 


— i.yo 


u.oo 


1 n 

iU 


— O.Uoi 


a 1 a 

o.io 


1 1 aa 
— i.ioo 


— i.yo 


u.oo 


ii 


— O.UOi 


a 
o.z 


1 1 aa 

— i.ioo 


— i.yo 


n 'iK 
u.oo 


1 

iz 


— O.i 


a naa 

o.Uoo 


1 1 aa 
— i.ioo 


— i.yo 


u.oo 




— O.Uo 


a naa 

o.Uoo 


1 1 aa 

— i.ioo 


— i.yo 


n Qc; 
u.oo 


1 /I 


— O.UO 


a naa 
o.Uoo 


1 1 aa 
— i. ioo 


1 QK 

— I.yo 


n Qc; 
u.oo 


iO 


0.U4I: 


a naa 

o.Uoo 


1 1 aa 

— 1. loo 


— i.yo 


u.oo 


1 

iO 


K no 
— o.uz 


a naa 
o.Uoo 


1 1 aa 
— 1. loo 


— i.yo 


u.oo 


1 7 


— 0. 


a naa 
o.Uoo 


1 1 aa 
— i.ioo 


— i.yo 


n "iK 
u.oo 


1 8 


— 4.yo 


a naa 

o.Uoo 


1 1 aa 

— i.ioo 


— i.yo 


n Qc; 
u.oo 


1 Q 


A QR 


a naa 
o.Uoo 


1 1 aa 
— i.ioo 


— i.yo 


u.oo 


on 
zu 


A OA 

— 4.y4 


a naa 

o.Uoo 


1 1 aa 
— i.ioo 


— i.yo 


u.oo 


21 


-4.92 


8.088 


-1.188 


-1.95 


0.35 


22 


-4.9 


8.088 


-1.188 


-1.95 


0.35 


23 


-5.031 


8.088 


-1.188 


-1.71 


0.35 


24 


-5.031 


8.088 


-1.188 


-1.76 


0.35 


25 


-5.031 


8.088 


-1.188 


-1.81 


0.35 


26 


-5.031 


8.088 


-1.188 


-1.86 


0.35 


27 


-5.031 


8.088 


-1.188 


-1.91 


0.35 


28 


-5.031 


8.088 


-1.188 


-1.96 


0.35 


29 


-5.031 


8.088 


-1.188 


-2. 


0.35 


30 


-5.031 


8.088 


-1.188 


-2.05 


0.35 



-2.3 



0.30 



0.32 



0.34 



0.36 



0.38 



0.40 



Fig. 1. — Relation between normalised polar moment of inertia c = C /niE? versus observed 
obliquity e in arcminutes. We compare the new for mula Eg .([21]) - black (thick) - with the 
usual one Eq.([T]) - red (dashed) - and with Eq.(4) of Peak (jl98ll ) - blue (dotted) - within in 
the interval 0.3 < c < 0.4 for the actual parameters of Mercury. 
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For each of them, a numerical study has been performed as explained above. After 
refinement of the initial conditions that gave us an accurate frequency decomposition of the 
signals, we identified 34 amplitudes, i.e. we can now write: 

K = J + fli — 2a2 cos {Q2 — ^1) + 2a3 cos (2^22 — 2Qi) — 20,4 cos {zui — VJ2) 

+ 2a5 cos {2wi — 2f2i) — 2a6 cos (3^2 — 3f2i) + 207 cos (4^2 — ^^i) 

+ 20,8 cos {w2 — wi + 0.2 — Oi) + 2a9 COS {wi — -072 + f22 — ^1) + 2aio cos {wi -\-W2 — 2O1) 

— 2aii cos {2wi — 30i + O2) — 2ai2 cos (5^2 — 50i) + 2ai3 cos {2wi — 2W2) 

— 2ai4 cos {vDi — W2 — 2Qi + 2^2) — 2ai5 cos {w2 — W\ — 2Jli + 2^2) 



(73 = 2ai7sin (Jl2 — ill) — 2ai8 sin (2Q2 — 2Qi) + 2ai9 sin (3Jl2 — 3Qi) 

— 2a2o sin (2roi — 2Qi) — 2a2i sin (4^2 — 4Qi) — 2022 sin {w2 — roi + Q2 — f^i) 

— 2a23 sin {w\ — ■072 + Jl2 — f^i) + 2a24 sin {^w\ — 3Qi + Q2) + 2a25 sin (5^2 — 5Qi) 

+ 2a26 sin (2roi — Qi — Q2) — 2027 sin {w\ + ZU2 — 2Qi) + 2a28 sin {—zui + ZU2 — 2Qi + 2^2) 

+ 2a29 sin {zui — ZU2 — 2Qi + 2Q2) — 2030 sin {2zui — 4Qi + 2^2) — 2031 sin {60,2 — 6O1) 

+ 2as2 sin {zui + ZU2 — 3Qi + O2) — 2033 sin {zui — ZU2 — 3Qi + 3O2) 

— 2034 sin {zu2 — w\ — 3Qi + 3^2) , (39) 



2ai6 cos (2^1 — 2Q2) , 



(38) 



with 



Cj (mR^) 



(40) 



' aC/ (mR^) + /3C20 + 7C22 + 5 
for i = 1, 2, 5, 6, 17, 18, 19, 20, 21, 22, 23, 24, 26, 27, 



C/ (mi?2) 



(41) 



q; + /3C20 + 7C22 



for i = 3,4 



C/ {mR^) 
a + I3C20 



(42) 



for i = 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 34, and 



C/ {mR^) 



(43) 



aC/ (mR^) + PC20 + 7 
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Table 3 — Continued 



1\T 


O20 (,X iU ) 




O30 (^XiU ) 


C-40 (,XiU j 


/ ymri ) 


31 


-5.031 


8.088 


-1.188 


-2.1 


0.35 


32 


-5.031 


8.088 


-1.188 


-2.15 


0.35 


33 


-5.031 


8.088 


-1.188 


-2.19 


0.35 


34 


-5.031 


8.088 


-1.188 


-1.95 


0.32 


35 


-5.031 


8.088 


-1.188 


-1.95 


0.326 


36 


-5.031 


8.088 


-1.188 


-1.95 


0.332 


37 


-5.031 


8.088 


-1.188 


-1.95 


0.338 


38 


-5.031 


8.088 


-1.188 


-1.95 


0.344 


39 


-5.031 


8.088 


-1.188 


-1.95 


0.349 


40 


-5.031 


8.088 


-1.188 


-1.95 


0.356 


41 


-5.031 


8.088 


-1.188 


-1.95 


0.362 


42 


-5.031 


8.088 


-1.188 


-1.95 


0.368 


43 


-5.031 


8.088 


-1.188 


-1.95 


0.374 


44 


-5.031 


8.088 


-1.188 


-1.95 


0.38 


45 


-5.031 


8.088 


-1.3 


-1.95 


0.35 


46 


-5.031 


8.088 


-1.28 


-1.95 


0.35 


47 


-5.031 


8.088 


-1.26 


-1.95 


0.35 


48 


-5.031 


8.088 


-1.24 


-1.95 


0.35 


49 


-5.031 


8.088 


-1.22 


-1.95 


0.35 


50 


-5.031 


8.088 


-1.2 


-1.95 


0.35 


51 


-5.031 


8.088 


-1.18 


-1.95 


0.35 


52 


-5.031 


8.088 


-1.16 


-1.95 


0.35 


53 


-5.031 


8.088 


-1.14 


-1.95 


0.35 


54 


-5.031 


8.088 


-1.12 


-1.95 


0.35 


55 


-5.031 


8.088 


-1.1 


-1.95 


0.35 
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for i = 25,28,29,30,31,32,33. 

The idea of these formulae come from Peale's formula (Eq|T]& Eq ET]) . We hoped to get 
amplitudes alike 



C/{mR^ 



' aCjimR^) + /3C20 + lC^2 + CC30 + <^C^^ + r/' 
For that, we tried to fit amplitudes with respect to one parameter, i.e. 



when possible and 



when not, and also 



a + hC- 



20 



a + hC- 



22 



(44) 



; I JL) ^ ^ "7 ('"^-). (45) 



/ ( ) = aC/ {mR') (46) 



/ {C20) = , (47) 



/(C22)= ,L ^ (4^ 



It turned out that it was impossible to estimate the influence of C30 and C40 with enough 
reliability, that is the reason why they do not appear in the final formulae (EqHOl to H3|) . 
The coefficients used are gathered in TabJH 



5. The influence of the different effects 



The derivation of these new formulae for the obliquity of Mercury allows us to estimate 
the influence of usually neglected effects, like the secular variations of the orbital elements. 
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the tides and the higher order harmonics. For that we first need to test the rehabihty of our 
initial conditions on a real, non-averaged simulation of the rotation of Mercury. 



5.1. Differences between the averaged and the unaveraged system 



We proceed as in (iNoyelles fc D'Hoedtl |2012| ). Sect. 4. To simulate the non-averaged 
rotation of Mercury, we integrate numerically the equations related to the following Hamil- 
tonian: 



77 



3GM 



ei (x' + y^) + €2 [x^ - y^) + eg z{bz^ - 3) + €4 (3 - 30^' + 35^^^ 



2 nd? 



with ei = -CiomR^/C, eg = 2C22mR^/C, €3 = C^quiR^ / {2C), 64 = C4omi?V(8C), 5 = 
1 — Cm/C, Cm being the polar inertial momentum of the mantle of Mercury, P is the norm 
of the angular momentum normalized by nC, d is the Sun-Mercury distance, and x and 
y are the first two coordinates of the unit vector pointing to the Sun in a reference frame 
defined by the principal axes of inertia of Mercury. The canonical variables associated with 
this Hamilltonian are 



/3, 



P 

R 



nCi 



the Hamilton equations associated being 



dh ^ &H 
dt aP' 

dh ^ dV. 
dt ai?' 



dP ^ _m 

dt dh ' 

dR ^ av. 

dt dh • 



The position of the Sun with respect to Mercury is computed using JPL DE406 ephemerides, 
so it contains every perturbation, i.e. not only the secular, but also the planetary ones. In- 
appropriate initial conditions in the numerical integration of the equations would yield un- 
expected free oscillations. The free longitudinal oscillations, their period being 12 years, 
can be easily damped adiabatically over 5,000 years, the DE406 ephemerides starting at 
-3,000. However, the free oscillations in obliquity, whose period is about 1,000 years, cannot 
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be damped over such a timescale without a significant and artificial impact on the equihb- 
rium. So, we add a damping only on the longitudinal motion, and we get free oscillations 
in obliquity as in Fig|2J This obliquity is the actual obliquity e, it is equal to K — I only if 
o"3 = 0. We obtain it with 



cos e 



G-n 



(52) 
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Fig. 2. — Obliquity e of Mercury given by the unaveraged system, with C; 



20 = -5.031 xlO'^ 

C22 = 8.088 X 10"^ C30 = C40 = and C = 0.35mR'^. The ^ 1,000 years oscillations are 
free librations that would not be present if the initial conditions were ideal. 

The amplitude of the free oscillations can be seen as an estimation of the error due 
to the initial conditions. This error can come from all the approximations made in the 
averaging process, in particular the limitation to a first order averaging, the expansions in 
eccentricity, or the exclusion of the planetary perturbations. The Figj3] shows the amplitude 
of these oscillations for different values of the interior parameters C20, C22 and C, the other 
ones not affecting our initial conditions. These amplitudes have been obtained thanks to a 
frequency analysis. W e can see that t he am plitude of the free oscillations are always smaller 
than 750 milli-arcsec. iMargot et al.l ( 12012| ) derived from observations a moment of inertia 
C = O.SAQmR^, so the amplitude of the free oscillations should be ^ 650 milli-arcsec. This 
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is the theoretical error induced by the Eql38] and Eqj39l 

niiiiiiiiiiiiiiiiiiiiiii 

C=0.39mR^ 

iiiiiiiKiiiiiiiiiniiiii 

C=0.35mR'^ 

llllillllllllllllllMIIII 

C=0.31mR^ 

I I I I I I I 

4 8.06 8.08 8.1 8.12 8.14 8.16 8.18 8.2 
C22 (x1e6) 

Fig. 3. — Amphtude of the free oscillations, induced by our initial conditions. 



5.2. Influence of J3 

We failed to find an influence of C30 = — J3 in our numerical and analytical formulae 
of the equilibrium obliquity. Anyway, this parameter is present in the full equations of the 
rotational dynamics of the system, and we checked its influence in plotting the average value 
of the obliquity e with respect to J3 (FigJl]). 

In FigJH we clearly see the influence of C30 on the obliquity e, which we did not see in 
the simple analytical nor the single averaged but still time dependent model. As we learned 
from (iii) of Sect. 2. 4, there is a time- dependent contribution proportional to C30 of the order 
of eR^/a^, which is affecting the action S3 - thus the inertial obliquity K as well as e through 
the presence of the resonant argument and the argument of the perihelion u. However, 
the effect is very small due to the presence of the in the denominator and can safely be 
neglected for the same reasoning we did to explain why J4 does not influence the results. 

A linear fit gives 



o 

w 
o 

1 

<n 
c 
o 

o 
w 
o 

<D 
03 



0.75 



0.7 - 



0.65 - 



0.6 



0.55 



mill 



iiiiiiii 



iiiiiiii 



7.98 8 8.02 8.0^ 



2.0658 



CO 



O 



o 



E 



2.0651 - 



2.0652 - 



2.0653 - 



2.0654 - 



2.0657 - 



2.0656 - 



2.0655 - 




2.065 



-1.3 -1.28 -1.26 -1.24 -1.22 -1.2 -1.18 -1.16 -1.14 -1.12 -1.1 

C30 (x1e5) 



Fig. 4. — Influence of C30 on tfie mean obliquity, obtained after numerical integration of the 
non-averaged equations. This has been computed from the cases 45 to 55 (Tabj3]). 



so neglecting J3 should induce an error of ~ 253 mas for J3 = 1.188 x 10 ^. We did the 
same job for J4 without finding any reliable influence. 

From a theoretical point of view the difference between the averaged (but still time 
dependent) and original system of equations of motion can be seen in the following way: 
with the average over the mean orbital longitude we neglect not only the fast periodic effects 
(the relevant timescale being the orbital period of Mercury), but we also set the orbital 
distance of Mercury from the Sun to a constant value. Thus, in the averaged system, we are 
unable to include all the perturbations, which affect the time evolution of the semimajor axis 
of the planet. In addition we needed to expand the potential into a truncated powerseries 
to be able to introduce the resonant argument and perform the averaging over the fast 
angle. Since we have limited all our expansions to the 4th order in eccentricity, an additional 
difference of the order C(e^) is expected. Last but not least, we use the simple average rule, 
which is equivalent to a first order averaging in the ratio of the masses m/M (the mass of 
Mercury over the mass of the Sun). 



e = 



(-355.I97C30 + 2.06115) arcmin. 



(53) 
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5.3. The tides 



The shape of Mercury, if hydrostatic, is due to the influence of a centrifugal potential, 
due to Mercury's spin, and a tidal potential. This tidal potential can be split into a static and 
an oscillating potential, as a consequence the gravity field parameters C-^n and C^? and the 
polar momentum of inertia C should experience periodic variations alike (IRappaport et al. 
20081 : 1 Van Hoolst et aljl2008l ): 



C20(t) 


/^static 1 
— '-^20 "T 


k2 
~2 


C22(t) 


/^static 
— ^22 ~ 


k2 
4 


m 


^static 


k2 
3 



(54) 
(55) 
(56) 



where ^2 is the classical Love number, e the eccentricity of Mercury, and U — I5 — Iq represents 
the mean anomaly of Mercury. We also have 



Qt 



m \ a 



(57) 



where M is the mass of the Sun, m the mass of Mercury, R its mean radius and a its 
semi-major axis. With e = 0.2056, a = 57,909,226.541 5 m (JPL HORIZONS), M/m = 
6.0239249 x 10^ and R = 2439.7 km dArchinal et aPboilh . we have g^e = -2.778369 x 10"^. 
The Love number ^2 should be between (fully inelastic rigid Mercury) and 1.5 (fully fluid 
Mercury). We think, from the detection of the longitudinal librations of Mercury, that it 
should be closer to a rigid body than to a fluid one, so we consider k2 = 0.5. The results 
are given in FigJSl We can see a peak-to-peak variation of ~ 100 mas. We can also see 
the secular drift due to the variations of the orbital elements related to eccentricity and 
inclination, ^ 10 mas over 20 years. 



5.4. Summary 



The influence of the different effects is gathered in the TabJS] The free librations that 
are mentioned are generated by the lack of accuracy of our initial conditions (Eqj38] fc [39|l . 
they do not have any physical relevance. We have made t he ass umptions that they have been 
damped to a negligible amplitude, as suggested by lPeald ( 120051 ). The other effects are smaller 
than that, they all have been estimated in this study except the polar motion, coming from 
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Fig. 5. — Variations of the obliquity of Mercury due to tides, with the numerical formula 
(left), and the analytical one (right). The period of variation is the orbital period, i.e. 88 
days. 



(INoyelles et al.l 120101 ). and the amplitude of the short-period librations, from (IDufey et al. 



20091 ). These numbers are very small compared to the accuracy of the determination of the 
obhquity, i.e. ~ 5 arcsec. 



6. Inverting the obliquity of Mercury 



Recently, iMargot et al.l (120121 ) measured an obliquity of (2.04 ± 0.08) arcmin, and they 
used Peale's formula (Eq{T]) to get C/{mR^) = 0.346 ± 0.014. In t his work, we presen t 
alternative formulae that we propose to use to invert the obliquity of iMargot et al.l (120121 ). 
We can consider that we propose 4 new formulae. From an analytical (Eq J^T]) and a numerical 
studies (Eqin21 from Eql3S] & EH]) we get 2 formulae. The influence of J3, that we detect 
only numerically with the full equations of the system (FigJH & Eq j53|) leads us to write 
down 2 new equations, in which the 2 obliquities given by the Eq|2T] & [52] are corrected by 
355.197 X J3. 

In considering that the spherical harmonics C20 = —5.031 x 10~^, C22 = 8.088 x 10^^, 
C30 = 1.188 X 10~^ and C40 = 1.95 x 10^^ are accurately known, we got the theoretical 
obliquity of Mercury 7 years after the date J2000.0 to be close to the mid-date of the radar 
observations. We considered that C was the only unknown parameter . In t he analytical 
formula, we used the same dynamical parameters as lYseboodt &: MargotI (120061 ) . i.e. t = 8.6° 
and a regressional period of the ascending node set to 328 kyr. And we get 



Margot et all (120121 ): C/{mR^) = 0.346 ± 0.014, 
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• Analytical formula (Eq^: C/{mR^) = 0.34621 ± 0.01358, 

• Numerical formula (Eqj52]): C/{mR^) = 0.34576 ± 0.01349, 

• Analytical formula (Eq|2l]) + J3: C/{mR^) = 0.34550 ± 0.01357, 

• Numerical formula (Eqj52]) + J3: C/{mR^) = 0.34506 ± 0.01348. 

All our numbers are consistent with the ones coming from Peale's formula. The polar 
moment of inertia of Mercury could be 0.345mi?^ instead of 0.346mi?^, this difference is 
ridiculous with respect to the accuracy of the observations. Giving so many digits lacks of 
physical relevance, but is necessary to stress the tiny differences between our 4 formulae. 



7. Conclusion 

This study tackles the influence of usually neglected effects like the secular variations 
of the orbital elements, the tides and the higher order harmonics on the instantaneous 
obliquity of Mercury. The main goal is to invert it to get clues on the internal structure, in 
particular the polar inertial momentum C. Moreover, it gives optimized initial conditions 
of the orientation of the angular momentum of Mercury, at any time and for any values of 
the internal struct ure parameters, that can b e directly used in numerical simulations. This 
is a refinement of (iNoyelles &: D'Hoedtll2012[ ). These initial conditions have the advantage 
to be Laplace Plane free, they are instead based on the ecliptic, whose definition is robust. 
They have been obtained thanks to averaged equations of the rotational motion, suitable 
for long-term studies. We hope that they will help fitting the r otation of Mercury by future 
experiments, like the radioscience experiment in BepiColombo ( iMilani et al.ll200ll ). 



We have shown that the usually neglected effects have an influence smaller than 1 arcsec, 
while the observations have an accuracy of ~ 5 arcsec. The determination of C can be at 
the most altered from 0.346mi?^ to 0.345mi?^, what does not fundamentally change our 
understanding of the internal structure of Mercury. The size of the molten core can be 
obtained from the long itudial librations of Me rcury, but depends on whether we consider a 
rigid inner core or not fIVan Hoolst et al.ll2012l ). 
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Appendix A 

We collect the various coefficients of the form [j]nm, with j,n,m G N, from Section [531 If 
we introduce the notation c^. = cos(x), = sin(x) they can be written in the form: 

At second order in (V20): 

[1]20 = 3 (1 + 3C2j) C2K , [2] 20 = ^SCiCRSiSK , [3] 20 = -GC2Ks'^ . 



This preprint was prepared with the AAS L^TJrjX macros v5.2. 
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At second order in {V22) 



[1]22 = 48424/2 , [2]22 = 12 (1 + q) (1 + Ck) S^Sk , [3]22 = 185,^4 , 
[4] 22 = ^^sl/^SiS^i^SK , [5] 22 = ^Ssf/2SK/2 , 



in (^20) and 



[6)22 = -44/2^? > [7)22 = 4 (1 + Ck) SiCiSK , [8)22 = - (1 + 3c2i) 4 > 
[9]22 = 4q {ck - 1) SiSK , [10)22 = -4s^4/2 



in (1*20) ■ At order three we have in ( V3 



30/ 



[IJao = -3042^' 4 , [2]30 = 30 {a - 1) (1 + 3ci) cks^s], , 
3 



[3]30 = -2 (13 + 20q + 15c2i) (3 + 5c2k) s\i^sk , 

3 3 
[4]3o = -| {^ck + 5c3k) {si + 5s3i) , [5]3o = (13 - 20q + 15c2i) {sk + 5S3k) , 

[6]3o = 30 (1 + a) {3ci - 1) CKSisl , [7]3o = -15 (q - 1) (1 + q) '4 . 
The fourth order coefficients, beeing part of (v^q), turn out to be: 

[1]40 = |r (9 + 20c2i + 35c4i) (9 + 20c2i^ + 35c4x) , [2^ = ^ (2s2i + 7s^i) {2s2k + 7s4k) , 
[3]4o = 15 (5 + 7c2i) (5 + 7c2k) s,'4 , [4]4o = 840QCKsf 4 , [5)40 = 105^4 . 

The fourth order coefficients in (^40) are: 

[6]4o = ¥ (5 + 7c2i) (9 + 20C2K + 35c4x) s^, , [7^0 = 840s1 8^4 , 
4 2 

[8]4o = 6720 (2q/2 + C3V2) cksI/^s% , 

[9]4o = 120 (9 + 14q + 7c2i) (5 + 7c2k) s-/24 , 

[10]40 = 30 (19Q/2 + 7 (2C3V2 + C5i/2)) sf/^ {2S2K + 7s^k) , 
[11]40 = -30^2 (l9Si/2 + 7 (s5i/2 - 2s3j/2)) (2s2X + 7s4k) , 

[12]4o = 120^2 (9 - 14q + 7c2i) (5 + 7c2k) 4 , 

[13]4o = 672042CK {s3i/2 - 2si/2) 4 , [14]4o = "210 {a - 1) (1 + Ci) ^4 . 
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Table 4: Coefficients involved in the EqHUlto H51 



N 


a 


/3 


7 


8 


1 


-9.6916394157 


-1.022791 X 10^ 


1.224118 X 10^ 


-0.041284174514 


2 


57.319667133 


-1.495053 X 10^ 


1.8067315 X 10^ 


-15.290923597 


3 


-288.588048 


-6.75794 X 10^ 


8.648745 x 10^ 


- 


4 


109.5839605 


-2.3327045 x 10^ 


-2.8315945 x 10^ 


- 


5 


9618.2490875 


-5.967955 x 10^ 


-1.4620515 X 10^° 


-4270.5575456 


6 


5690.3083952 


-3.29021 X 10^ 


4.59872 X 10^ 


-5791.2841683 


7 


160882.75 


-1.563842 X 10^° 


- 


- 


8 


-330146.25 


-3.3956545 x 10^° 


- 


- 


9 


-307545.35 


-3.441592 X 10^° 


- 


- 


10 


-879441.5 


-4.59746 X 10^° 


- 


- 


11 


-1025209.5 


-5.12722 X 10^° 


- 


- 


12 


760424 


-7.318535 X 10^° 


- 


- 


13 


-3583265 


-1.7097535 X 10^^ 


- 


- 


14 


-1352596 


-1.522164 X 10^1 


- 


- 


15 


-1475953.5 


-1.554497 X 10^1 


- 


- 


16 


-4788455 


-2.4530905 x 10^^ 


- 


- 


17 


0.0459703076 


-111936.3 


133840.35 


0.00068014043987 


18 


0.5097512707 


-474698 


568708 


-0.0063487867442 


19 


3.6758306831 


-2005542 


2431534 


-0.32882074509 


20 


20.470309592 


-12322835 


-31382330 


1.6246447377 


21 


22.351343806 


-8470000 


10602865 


-4.3502813922 


22 


1.8968370715 


-25710160 


-28807205 


-14.991580755 


23 


28.007891612 


-25704595 


-30518425 


0.75276141087 


24 


151.69560968 


-52410050 


-127924300 


-54.254151937 


25 


11.320467298 


-35729050 


384.87733645 


- 


26 


60.949766585 


-91723450 


-262171000 


275.67522095 


27 


93.149257619 


-96422550 


-265414800 


257.09346578 


28 


256.01109782 


-108893750 


-1098.8078842 




29 


93.896968203 


-109386550 


-1017.8654389 




30 


880.39601125 


-223772150 


-4722.8636039 




31 


-390.02956574 


-145710250 


2105.3653480 




32 


1929.3930214 


-404092500 


-8327.8625575 




33 


2613.8902647 


-456302000 


-4855.3715926 




34 


-2905.1715 


-438245500 
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Table 5: Relative influence on the mean obliquity of usually neglected effects. The amplitude 
of the free librations can be seen as an estimation of the error due to theory. 



Effect 


Influence on obhquity 


Free librations 


< 750 mas 


C30 


^ 250 mas 


Tides 


^ 100 mas 


Polar motion 


^ 80 mas fNovelles et al. 20101 


Short-period librations 


< 20 mas fDufev et al. 2009) 


Secular drift 


^ 10 mas over 20 years 


C40 


negligible 



